CRUK CI Summer School 2021

Pseudotime Analysis

1 Pseudotime analysis

1.1 Learning objectives

library(SingleCellExperiment)
library(scran)
library(scater)
library(batchelor)
library(cowplot)
library(pheatmap)
library(tidyverse)
library(SingleR)
library(destiny)
library(gam)
library(viridis)
library(msigdbr)
library(clusterProfiler)
#library(cellAlign) # https://github.com/shenorrLab/cellAlign
#library(Cairo)

1.2 Extract T-cells from HCA ABMMC Dataset

In this section, we are starting our analysis with normalized HCA data and perform integration, clustering and dimensionality reduction. Our aim is to extract T-cells from this dataset and proceed with pseudotime analysis in the next section.

We are going to work with HCA data. This data set has been pre-processed and normalized before.

sce <- readRDS(\../CourseMaterials/Robjects/hca_sce_nz_postDeconv_5kCellPerSpl.Rds\)

We use symbols in place of ENSEMBL IDs for easier interpretation later.

rownames(sce) <- uniquifyFeatureNames(rowData(sce)$ensembl_gene_id, names = rowData(sce)$Symbol)

1.2.1 Variance modeling

We block on the donor of origin to mitigate batch effects during highly variable gene (HVG) selection. We select a larger number of HVGs to capture any batch-specific variation that might be present.

dec.hca <- modelGeneVar(sce, block=sce$Sample.Name)
top.hca <- getTopHVGs(dec.hca, n=5000)

1.2.2 Data integration

The batchelor package provides an implementation of the Mutual Nearest Neighbours (MNN) approach via the fastMNN() function. We apply it to our HCA data to remove the donor specific effects across the highly variable genes in top.hca. To reduce computational work and technical noise, all cells in all samples are projected into the low-dimensional space defined by the top d principal components. Identification of MNNs and calculation of correction vectors are then performed in this low-dimensional space. The corrected matrix in the reducedDims() contains the low-dimensional corrected coordinates for all cells, which we will use in place of the PCs in our downstream analyses. We store it in ‘MNN’ slot in the main sce object.

set.seed(1010001)
merged.hca <- fastMNN(sce,
              batch = sce$Sample.Name,
              subset.row = top.hca)
reducedDim(sce, 'MNN') <- reducedDim(merged.hca, 'corrected')

1.2.3 Dimensionality Reduction

We cluster on the low-dimensional corrected coordinates to obtain a partitioning of the cells that serves as a proxy for the population structure. If the batch effect is successfully corrected, clusters corresponding to shared cell types or states should contain cells from multiple samples. We see that all clusters contain contributions from each sample after correction.

set.seed(01010100)
sce <- runPCA(sce, dimred=\MNN\)
sce <- runUMAP(sce, dimred=\MNN\)
sce <- runTSNE(sce, dimred=\MNN\)

plotPCA(sce, colour_by=\Sample.Name\) + ggtitle(\PCA\)

plotTSNE(sce, colour_by=\Sample.Name\) + ggtitle(\tSNE\)

plotUMAP(sce, colour_by=\Sample.Name\) + ggtitle(\UMAP\)

1.2.4 Clustering

Graph-based clustering generates an excessively large intermediate graph so we will instead use a two-step approach with k-means. We generate 1000 small clusters that are subsequently aggregated into more interpretable groups with a graph-based method.

set.seed(1000)
clust.hca <- clusterSNNGraph(sce,
                             use.dimred=\MNN\,
                             use.kmeans=TRUE,
                             kmeans.centers=1000)
colLabels(sce) <- factor(clust.hca)
table(colLabels(sce))

    1     2     3     4     5     6     7     8     9    10    11 
 4471  6365  2911  4420  2402 12536  1176  3566  1146   596   411 

plotPCA(sce, colour_by=\label\) + ggtitle(\PCA\)

plotUMAP(sce, colour_by=\label\) + ggtitle(\UMAP\)

plotTSNE(sce, colour_by=\label\) + ggtitle(\tSNE\)

1.2.5 Cell type classification

We perform automated cell type classification using a reference dataset to annotate each cluster based on its pseudo-bulk profile. This is for a quick assignment of cluster identity. We are going to use Human Primary Cell Atlas (HPCA) data for that. HumanPrimaryCellAtlasData function provides normalized expression values for 713 microarray samples from HPCA (Mabbott et al., 2013). These 713 samples were processed and normalized as described in Aran, Looney and Liu et al. (2019). Each sample has been assigned to one of 37 main cell types and 157 subtypes.

se.aggregated <- sumCountsAcrossCells(sce, id=colLabels(sce))
hpc <- celldex::HumanPrimaryCellAtlasData()
anno.hca <- SingleR(se.aggregated, ref = hpc, labels = hpc$label.main, assay.type.test=\sum\)
anno.hca

DataFrame with 11 rows and 5 columns
                           scores     first.labels       tuning.scores
                         <matrix>      <character>         <DataFrame>
1  0.270340:0.754119:0.611547:...           B_cell 0.754119: 0.6997738
2  0.275371:0.598736:0.693408:... Pre-B_cell_CD34- 0.509690: 0.0634099
3  0.364176:0.613523:0.692287:...              CMP 0.594369: 0.3957973
4  0.282579:0.620968:0.575578:...          NK_cell 0.584395: 0.4455476
5  0.381883:0.558646:0.649733:...              MEP 0.361367: 0.3372588
6  0.291896:0.637409:0.575667:...          T_cells 0.773909: 0.7208437
7  0.272533:0.650300:0.600539:...          NK_cell 0.805271: 0.7248806
8  0.296865:0.637337:0.590623:...          T_cells 0.689381:-0.0530141
9  0.315724:0.599837:0.722784:... Pre-B_cell_CD34- 0.475323: 0.3332805
10 0.321822:0.686603:0.592921:...           B_cell 0.488044: 0.2851121
11 0.293917:0.635528:0.597594:... Pre-B_cell_CD34- 0.245206: 0.2048503
             labels    pruned.labels
        <character>      <character>
1            B_cell           B_cell
2          Monocyte         Monocyte
3               CMP              CMP
4           T_cells          T_cells
5        BM & Prog.       BM & Prog.
6           T_cells          T_cells
7           NK_cell          NK_cell
8           T_cells          T_cells
9  Pre-B_cell_CD34- Pre-B_cell_CD34-
10           B_cell           B_cell
11 Pre-B_cell_CD34- Pre-B_cell_CD34-

tab <- table(anno.hca$labels, colnames(se.aggregated))
# Adding a pseudo-count of 10 to avoid strong color jumps with just 1 cell.
pheatmap(log10(tab+10))

sce$cell_type<-recode(sce$label,
       "1" = "T_cells", 
       "2" = "Monocyte", 
       "3"="B_cell",
       "4"="MEP", 
       "5"="B_cell", 
       "6"="CMP", 
       "7"="T_cells",
      "8"="Monocyte",
      "9"="T_cells",
      "10"="Pro-B_cell_CD34+",
      "11"="NK_cell",
      "12"="B_cell")

#level_key <- anno.hca %>%
#  data.frame() %>%
#  rownames_to_column(\clu\) %>%
#  #select(clu, labels)
#  pull(labels)
level_key <- anno.hca$labels
names(level_key) <- row.names(anno.hca)
sce$cell_type <- recode(sce$label, !!!level_key)

We can now use the predicted cell types to color PCA, UMAP and tSNE.

plotPCA(sce, colour_by=\cell_type\, text_by=\cell_type\) + ggtitle(\PCA\)

plotUMAP(sce, colour_by=\cell_type\, text_by=\cell_type\) + ggtitle(\UMAP\)

plotTSNE(sce, colour_by=\cell_type\, text_by=\cell_type\) + ggtitle(\tSNE\)

We can also check expression of some marker genes.

CD3D and TRAC are used as marker genes for T-cells Szabo et al. 2019.

plotExpression(sce, features=c(\CD3D\), x=\label\, colour_by=\cell_type\)

plotExpression(sce, features=c(\TRAC\), x=\label\, colour_by=\cell_type\)

1.2.6 Extract T-cells

We will now extract T-cells and store in a new SCE object to use in pseudotime analysis.

Pull barcodes for T-cells

tcell.bc <- colData(sce) %>%
    data.frame() %>%
    group_by(cell_type) %>%
    dplyr::filter(cell_type == \T_cells\) %>%
    pull(Barcode)

table(colData(sce)$Barcode %in% tcell.bc)

FALSE  TRUE 
19478 20522 

Create a new SingleCellExperiment object for T-cells

tmpInd <- which(colData(sce)$Barcode %in% tcell.bc)
sce.tcell <- sce[,tmpInd]

saveRDS(sce.tcell,\../CourseMaterials/Robjects/sce.tcell.RDS\)
#rm(sce.tcell)

1.3 Setting up the data

In many situations, one is studying a process where cells change continuously. This includes, for example, many differentiation processes taking place during development: following a stimulus, cells will change from one cell-type to another. Ideally, we would like to monitor the expression levels of an individual cell over time. Unfortunately, such monitoring is not possible with scRNA-seq since the cell is lysed (destroyed) when the RNA is extracted.

Instead, we must sample at multiple time-points and obtain snapshots of the gene expression profiles. Since some of the cells will proceed faster along the differentiation than others, each snapshot may contain cells at varying points along the developmental progression. We use statistical methods to order the cells along one or more trajectories which represent the underlying developmental trajectories, this ordering is referred to as “pseudotime”.

A recent benchmarking paper by Saelens et al provides a detailed summary of the various computational methods for trajectory inference from single-cell transcriptomics. They discuss 45 tools and evaluate them across various aspects including accuracy, scalability, and usability. They provide dynverse, an open set of packages to benchmark, construct and interpret single-cell trajectories (currently they have a uniform interface for 60 methods). We load the SCE object we have generated previously. This object contains only the T-cells from 8 healthy donors. We will first prepare the data by identifying variable genes, integrating the data across donors and calculating principal components.

sce.tcell <- readRDS(\../CourseMaterials/Robjects/sce.tcell.RDS\)
sce.tcell

class: SingleCellExperiment 
dim: 20425 20522 
metadata(0):
assays(2): counts logcounts
rownames(20425): AL627309.1 AL669831.5 ... AC233755.1 AC240274.1
rowData names(11): ensembl_gene_id external_gene_name ... detected
  gene_sparsity
colnames: NULL
colData names(19): Sample Barcode ... label cell_type
reducedDimNames(4): MNN PCA UMAP TSNE
altExpNames(0):

dec.tcell <- modelGeneVar(sce.tcell, block=sce.tcell$Sample.Name)
top.tcell <- getTopHVGs(dec.tcell, n=5000)

set.seed(1010001)
merged.tcell <- fastMNN(sce.tcell, batch = sce.tcell$Sample.Name, subset.row = top.tcell)
reducedDim(sce.tcell, 'MNN') <- reducedDim(merged.tcell, 'corrected')

sce.tcell <- runPCA(sce.tcell, dimred=\MNN\)

plotPCA(sce.tcell, colour_by=\Sample.Name\)

1.3.1 Trajectory inference with destiny

Diffusion maps were introduced by Ronald Coifman and Stephane Lafon, and the underlying idea is to assume that the data are samples from a diffusion process. The method infers the low-dimensional manifold by estimating the eigenvalues and eigenvectors for the diffusion operator related to the data. Angerer et al have applied the diffusion maps concept to the analysis of single-cell RNA-seq data to create an R package called destiny.

For ease of computation, we will perform pseudotime analysis only on one sample, and we will downsample the object to 1000 cells. We will select the sample named MantonBM1.

# pull the barcodes for MantonBM1 sample & and downsample the set to 1000 genes 
vec.bc <- colData(sce.tcell) %>%
    data.frame() %>%
    filter(Sample.Name == \MantonBM1\) %>%
    group_by(Sample.Name) %>%
    sample_n(1000) %>%
    pull(Barcode)

Number of cells in the sample:

table(colData(sce.tcell)$Barcode %in% vec.bc)

FALSE  TRUE 
19522  1000 

Subset cells from the main SCE object:

tmpInd <- which(colData(sce.tcell)$Barcode %in% vec.bc)
sce.tcell.BM1 <- sce.tcell[,tmpInd]
sce.tcell.BM1

class: SingleCellExperiment 
dim: 20425 1000 
metadata(0):
assays(2): counts logcounts
rownames(20425): AL627309.1 AL669831.5 ... AC233755.1 AC240274.1
rowData names(11): ensembl_gene_id external_gene_name ... detected
  gene_sparsity
colnames: NULL
colData names(19): Sample Barcode ... label cell_type
reducedDimNames(4): MNN PCA UMAP TSNE
altExpNames(0):

Identify top 500 highly variable genes

dec.tcell.BM1 <- modelGeneVar(sce.tcell.BM1)
top.tcell.BM1 <- getTopHVGs(dec.tcell.BM1, n=500)

We will extract normalized counts for HVG to use in pseudotime alignment

tcell.BM1_counts <- logcounts(sce.tcell.BM1)
tcell.BM1_counts <- t(as.matrix(tcell.BM1_counts[top.tcell.BM1,]))
cellLabels <- sce.tcell.BM1$Barcode
rownames(tcell.BM1_counts) <- cellLabels

tcell.BM1_counts[1:4,1:4]

                        CCL5     NKG7   S100A4 KLRB1
AAACCTGAGCAGCGTA-13 3.733773 1.609114 2.838333     0
AAACCTGAGCGATATA-13 0.000000 0.000000 3.889855     0
AAACCTGGTAATCACC-13 0.000000 0.000000 3.576974     0
AAACCTGGTACAGTTC-13 3.337079 1.712109 3.630169     0

And finally, we can run pseudotime alignment with destiny

dm_tcell_BM1 <- DiffusionMap(tcell.BM1_counts,n_pcs = 50)

Plot diffusion component 1 vs diffusion component 2 (DC1 vs DC2).

tmp <- data.frame(DC1 = eigenvectors(dm_tcell_BM1)[, 1],
                  DC2 = eigenvectors(dm_tcell_BM1)[, 2])

ggplot(tmp, aes(x = DC1, y = DC2)) +
    geom_point() + 
    xlab(\Diffusion component 1\) + 
    ylab(\Diffusion component 2\) +
    theme_classic()

Stash diffusion components to SCE object

sce.tcell.BM1$pseudotime_destiny_1 <- eigenvectors(dm_tcell_BM1)[, 1]
sce.tcell.BM1$pseudotime_destiny_2 <- eigenvectors(dm_tcell_BM1)[, 2]

1.3.2 Find temporally expressed genes

After running destiny, an interesting next step may be to find genes that change their expression over the course of time We demonstrate one possible method for this type of analysis on the 500 most variable genes. We will regress each gene on the pseudotime variable we have generated, using a general additive model (GAM). This allows us to detect non-linear patterns in gene expression. We are going to use HVG we identified in the previous step, but this analysis can also be done using the whole transcriptome.

# Only look at the 500 most variable genes when identifying temporally expressesd genes.
# Identify the variable genes by ranking all genes by their variance.
# We will use the first diffusion components as a measure of pseudotime 
Y<-log2(counts(sce.tcell.BM1)+1)
colnames(Y)<-cellLabels
Y<-Y[top.tcell.BM1,]
# Fit GAM for each gene using pseudotime as independent variable.
t <- eigenvectors(dm_tcell_BM1)[, 1]
gam.pval <- apply(Y, 1, function(z){
  d <- data.frame(z=z, t=t)
  tmp <- gam(z ~ lo(t), data=d)
  p <- summary(tmp)[4][[1]][1,5]
  p
})

Select top 30 genes for visualization

# Identify genes with the most significant time-dependent model fit.
topgenes <- names(sort(gam.pval, decreasing = FALSE))[1:30]  

Visualize these genes in a heatmap

heatmapdata <- Y[topgenes,]
heatmapdata <- heatmapdata[,order(t, na.last = NA)]
t_ann<-as.data.frame(t)
colnames(t_ann)<-\pseudotime\
pheatmap(heatmapdata, cluster_rows = T, cluster_cols = F, color = plasma(200), show_colnames = F, annotation_col = t_ann)

Visualize how some of the temporally expressed genes change in time

Following individual genes is very helpful for identifying genes that play an important role in the differentiation process. We illustrate the procedure using the GZMA gene. We have added the pseudotime values computed with destiny to the colData slot of the SCE object. Having done that, the full plotting capabilities of the scater package can be used to investigate relationships between gene expression, cell populations and pseudotime.

plotExpression(sce.tcell.BM1,
           \GZMA\,
           x = \pseudotime_destiny_1\, 
               show_violin = TRUE,
               show_smooth = TRUE)

1.3.3 Pseudotime analysis for another HCA sample

# pull barcodes for MantonBM2 
vec.bc <- colData(sce.tcell) %>%
    data.frame() %>%
    filter(Sample.Name == "MantonBM2") %>%
    group_by(Sample.Name) %>%
    sample_n(1000) %>%
    pull(Barcode)

Error in colData(sce.tcell) : object 'sce.tcell' not found